---
title: "Overcoming Impeachment Hurdles: Elite Polarization, Mass Mobilization and Corruption Scandals"
paper_authors: "Mahmoud Farag, Isabella C. Montini and Philipp Schemm"
script_author: "Philipp Schemm"
date: "2026-03-03"
output:
  pdf_document:
    toc: yes
    toc_depth: '3'
  html_document:
    df_print: paged
    toc: yes
    toc_depth: '3'
---

# Set directory
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
getwd()
path <- setwd() #insert your directory
```


# Data preparation
```{r}
# load packages 
library(QCA)
library(SetMethods)
library(readxl)
library(readr)
library(writexl)
library(data.table)

list.files()
data_impeach <- read.csv("data_impeach_main.csv")
colnames(data_impeach)

# raw data 
attach(data_impeach)  
data_r <- data.frame(code,
                     out_lower_chamber_vote,
                     out_removal, 
                     mass_mobilization, 
                     leg_shield, 
                     elite_ide_polar, 
                     elite_aff_polar, 
                     corrupt_scandal, 
                     e_regionpol_7C) 
detach(data_impeach)
```

# Calibration of set data (see full set data in Appendix 3)
```{r}
# In the paper, we further shortened the variable names to ensure readability.
# Caution: In the paper, we refer to lower chamber vote as outcome 1 and removal from office as outcome 2. However, in this script, some objects include the name out2 and out3, which refer respectively to out_lower_chamber_vote and out_removal. For example, "ttout2" is the truth table for out_lower_chamber_vote, while "ttout3" is the truth table for the outcome out_removal. 

# analysis                  paper
# out_lower_chamber_vote  = out1, 
# out_removal             = out2, 
# elite_ide_polar         = elip, 
# elite_aff_polar         = elap, 
# leg_shield              = les, 
# corrupt_scandal         = csca, 
# mass_mobilization       = mam

df <- data.frame(
  code = data_r$code,                                       # case-identificator
  e_regionpol_7C = data_r$e_regionpol_7C,                   # regional identificator for cluster analysis
  out_lower_chamber_vote = data_r$out_lower_chamber_vote,   # out1
  out_removal = data_r$out_removal,                         # out2
  elite_ide_polar = calibrate(data_r$elite_ide_polar,       # elip
                              thresholds = 5, 
                              type = "crisp"),
  elite_aff_polar = calibrate(4-data_r$elite_aff_polar,     # elap  ## "4-" because it is a reversed scale: This way "always", "usual" and "about half of the time" are coded 1, "usually not" and  "never" are coded 0. Cases with exactly 1 (!3 in the non-reversed raw data!) should be coded as 0, so we have to use a value sightly higher than 1.
                              thresholds = 1.001,           
                              type = "crisp"),
  mass_mobilization = calibrate(data_r$mass_mobilization,   # mam
                                thresholds = 10000, 
                                type = "crisp"),
  leg_shield = calibrate(data_r$leg_shield,                 # les
                         thresholds = 0.1, 
                         type = "crisp"),
  corrupt_scandal = data_r$corrupt_scandal                  # csca
)

# subset set data for the two outcomes
df_set_out1 <- df 
df_set_out2 <- subset(df, out_lower_chamber_vote == 1)

# assign president-years as rownames
df_set_out1 <- data.frame(df_set_out1[,-1], row.names=df_set_out1[,1])
df_set_out2 <- data.frame(df_set_out2[,-1], row.names=df_set_out2[,1])


```

# Skewness tests for set data (Appendix 4)
```{r}
# function for skewness checks
check_skewness <- function(data) {
  # extract numeric columns
  numeric_cols <- sapply(data, is.numeric)
  if (sum(numeric_cols) == 0) {
    message("Keine numerischen Spalten im Data Frame gefunden.")
    return(NULL)
  }
  
  # step 2: calculate skewness
  skew_vector <- sapply(data[, numeric_cols], SetMethods::skew.check)
  
  # step 3: creating a table
  skew_df <- data.frame(
    Variable = names(skew_vector),
    Schiefe = skew_vector,
    row.names = NULL  
  )
  
  return(skew_df)
}

# select only relevant variables
df_set_out1_skw <- subset(df_set_out1, select = c(out_lower_chamber_vote, elite_ide_polar, elite_aff_polar, leg_shield, corrupt_scandal, mass_mobilization))

df_set_out2_skw <- subset(df_set_out2, select = c(out_removal, elite_ide_polar, elite_aff_polar, leg_shield, corrupt_scandal, mass_mobilization))
  
# check skewness
check_skewness(df_set_out1_skw) # out_lower_chamber_vote
check_skewness(df_set_out2_skw) # out_removal

```

# Necessity analysis (Appendix 5 and Appendix 9)
```{r}

# Outcome: out_lower_chamber_vote

# check for single necessary conditions

cons_out2 <- subset(df_set_out1, select = c("elite_ide_polar", "elite_aff_polar","mass_mobilization", "leg_shield", "corrupt_scandal") )

pof(cons_out2, "out_lower_chamber_vote", df_set_out1, relation = "nec") # presence of conditions and outcome
pof(1-cons_out2, "out_lower_chamber_vote", df_set_out1, relation = "nec") # absence of conditions
pof(cons_out2, "~out_lower_chamber_vote", df_set_out1, relation = "nec") # absence of outcome
pof(1-cons_out2, "~out_lower_chamber_vote", df_set_out1, relation = "nec") # absence of conditions and outcome
#no single necessary condition could be found, neither for the outcome, nor for the non-outcome

#check for SUIN condition
# SUIN_out2 <- superSubset(df_set_out1, outcome = "out_lower_chamber_vote", conditions = "elite_ide_polar, elite_aff_polar, mass_mobilization, leg_shield, corrupt_scandal", relation = "nec", incl.cut = 0.9, ron.cut = 0.5, cov.cut = 0.6)
# SUIN_out2 # no SUIN condition could be found for the outcome

#absence of outcome
SUIN_noout2 <- superSubset(df_set_out1, outcome = "~out_lower_chamber_vote", conditions = "elite_ide_polar, elite_aff_polar, mass_mobilization, leg_shield, corrupt_scandal", relation = "nec", incl.cut = 0.9, ron.cut = 0.5, cov.cut = 0.6)
SUIN_noout2 # no SUIN condition could be found for the non-outcome


# Outcome: out_removal

# check for single necessary conditions

cons_out3 <- subset(df_set_out2, select = c("elite_ide_polar", "elite_aff_polar","mass_mobilization", "leg_shield", "corrupt_scandal") )

pof(cons_out3, "out_removal", df_set_out2, relation = "nec")
pof(1-cons_out3, "out_removal", df_set_out2, relation = "nec")
pof(cons_out3, "~out_removal", df_set_out2, relation = "nec")
pof(1-cons_out3, "~out_removal", df_set_out2, relation = "nec")
#no single necessary condition could be found, neither for the outcome, nor for the non-outcome

#check for SUIN condition
SUIN_out3 <- superSubset(df_set_out2, outcome = "out_removal", conditions = "elite_ide_polar, elite_aff_polar, mass_mobilization, leg_shield, corrupt_scandal", relation = "nec", incl.cut = 0.9, ron.cut = 0.5, cov.cut = 0.6)
SUIN_out3 # one candidate SUIN condition could be found for the outcome, but it has low RoN (0.5)  

#absence of outcome
SUIN_noout3 <- superSubset(df_set_out2, outcome = "~out_removal", conditions = "elite_ide_polar, elite_aff_polar, mass_mobilization, leg_shield, corrupt_scandal", relation = "nec", incl.cut = 0.9, ron.cut = 0.5, cov.cut = 0.6)
SUIN_noout3 # no SUIN condition could be found for the outcome

```

# Sufficiency Analysis
## Truth table 
```{r}

# truth tables for out_lower_chamber_vote
ttout2 <- truthTable(data=df_set_out1, outcome = "out_lower_chamber_vote", conditions = "elite_ide_polar, elite_aff_polar, mass_mobilization, leg_shield, corrupt_scandal", 
                          incl.cut=0.9,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE, 
                          dcc = TRUE)
ttout2 

# truth tables for out_removal 
ttout3 <- truthTable(data=df_set_out2, outcome = "out_removal", conditions = "elite_ide_polar, elite_aff_polar, mass_mobilization, leg_shield, corrupt_scandal", 
                          incl.cut=1,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE,
                          dcc = TRUE)
ttout3 

```
## Minimization (conservative, parsimonious and intermediate solution)
```{r}

## Outcome: out_lower_chamber_vote

# conservative solution csttout2
csttout2 <- minimize(ttout2, details=TRUE, show.cases=TRUE, row.dom=TRUE, all.sol=TRUE, use.tilde=TRUE) 
csttout2
csttout2$PIchart

# parsimonious solution psttout2 
psttout2 <- minimize(ttout2, include="?", details=TRUE, show.cases=TRUE, row.dom=TRUE, all.sol=TRUE, use.tilde=TRUE)
psttout2
psttout2$PIchart 
psttout2$SA

# intermediate solution isttout2 (with directional expectations)
isttout2 <- minimize(ttout2, include="?", dir.exp="-,-,1,0,1", details=TRUE, show.cases=TRUE, row.dom=TRUE,
                         all.sol=TRUE, use.tilde=TRUE)
isttout2
isttout2$PIchart
isttout2$i.sol$C1P1$EC 
isttout2$i.sol$C1P1$DC 

## outcome: out_removal

# conservative solution csttout3
csttout3 <- minimize(ttout3, details=TRUE, show.cases=TRUE, row.dom=TRUE, all.sol=TRUE, use.tilde=TRUE) 
csttout3
csttout3$PIchart #prime implicants chart

# parsimonious solution psttout3 
psttout3 <- minimize(ttout3, include="?", details=TRUE, show.cases=TRUE, row.dom=TRUE, all.sol=TRUE, use.tilde=TRUE)
psttout3
psttout3$PIchart 
psttout3$SA

# intermediate solution isttout3 (with directional expectations)
isttout3 <- minimize(ttout3, include="?", dir.exp="-,-,1,0,1", details=TRUE, show.cases=TRUE, row.dom=TRUE,
                         all.sol=TRUE, use.tilde=TRUE)
isttout3
isttout3$PIchart
isttout3$i.sol$C1P1$EC #easy counterfactuals
isttout3$i.sol$C1P1$DC #difficult counterfactuals

```

# Enhanced sufficiency analysis for the outcome (includes results shown in Table 1-4)
```{r}

## Outcome: out_lower_chamber_vote 

# Steps:
#1 Contradictions to necessity:
# no necessary conditions for the outcome -> no contradictions possible
#2 Implausible assumptions
# There is no possible implausible assumption in my opinion. 
#3 Simultaneous subset relations 
ttnoout2 <- truthTable(data=df_set_out1, outcome = "~out_lower_chamber_vote", conditions = "elite_ide_polar, elite_aff_polar, mass_mobilization, leg_shield, corrupt_scandal", 
                          incl.cut=0.9,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE, 
                          dcc = TRUE)
ttnoout2
# compare if any of the rows that you included into the solution is also part of a possible solution for the non-outcome
ttnoout2
ttout2

#either look manually. No simultaneos subset relations found. Or with formula:
SSR <- intersect(rownames(ttout2$tt) [ttout2$tt$OUT==1],
                                rownames(ttnoout2$tt) [ttnoout2$tt$OUT==1])
SSR
# = 0
# -> no simultaneous subset relations can be found
#4 Incoherent counterfactuals
contcountercrisp2_out2 <- findRows(obj=ttout2, type = 2)
contcountercrisp2_out2
# Incoherent rows: 4 12 16 20 27

#to be able to compare them, a new truth table is constructed
ttout2enhanced <- truthTable(data=df_set_out1, outcome = "out_lower_chamber_vote", conditions = "elite_ide_polar, elite_aff_polar, mass_mobilization, leg_shield, corrupt_scandal", 
                          incl.cut=0.9,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE, 
                          dcc = TRUE)
ttout2enhanced

# turn remainder rows to 0
ttout2enhanced$tt['4', 'OUT'] <- 0
ttout2enhanced$tt['12', 'OUT'] <- 0
ttout2enhanced$tt['16', 'OUT'] <- 0
ttout2enhanced$tt['20', 'OUT'] <- 0
ttout2enhanced$tt['27', 'OUT'] <- 0

#enhanced truth table
ttout2enhanced
#write.csv(ttout2enhanced$tt, "ttout2enhanced.csv")

#enhanced parsimonious solution
psttout2enhanced <- minimize(ttout2enhanced, include="?", details=TRUE, show.cases=TRUE, row.dom=TRUE, all.sol=TRUE, use.tilde=TRUE)
psttout2enhanced 

psttout2 # for comparison

#enhanced intermediate solution
isttout2enhanced <- minimize(ttout2enhanced, include="?", dir.exp="-,-,1,0,1", details=TRUE, show.cases=TRUE, row.dom=TRUE,
                         all.sol=TRUE, use.tilde=TRUE)
isttout2enhanced  ######### FINAL SOLUTION FOR THE OUTCOME "lower house vote" #########
isttout2enhanced$i.sol$C1P1$EC 
isttout2enhanced$i.sol$C1P1$DC

isttout2 # original version for comparison



## Outcome: out_removal

# Steps:
#1 Contradictions to necessity:
# no necessary conditions for the outcome -> no contradictions possible
#2 Implausible assumptions
# There is no possible implausible assumption. 
#3 Simultaneous subset relations 
ttnoout3 <- truthTable(data=df_set_out2, outcome = "~out_removal", conditions = "elite_ide_polar, elite_aff_polar, mass_mobilization, leg_shield, corrupt_scandal", 
                          incl.cut=1,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE,
                          dcc = FALSE)
ttnoout3
# compare if any of the rows that you included into the solution is also part of a possible solution for the non-outcome
ttnoout3
ttout3

#either look manually. No simultaneos subset relations found. Or with formula:
SSR <- intersect(rownames(ttout3$tt) [ttout3$tt$OUT==1],
                                rownames(ttnoout3$tt) [ttnoout3$tt$OUT==1])
SSR
# = 0
# -> no simultaneous subset relations can be found
#4 Incoherent counterfactuals
contcountercrisp2_out3 <- findRows(obj=ttout3, type = 2)
contcountercrisp2_out3
# Incoherent rows: 3  7 19 23

#to be able to compare them, a new truth table is constructed
ttout3enhanced <- truthTable(data=df_set_out2, outcome = "out_removal", conditions = "elite_ide_polar, elite_aff_polar, mass_mobilization, leg_shield, corrupt_scandal", 
                          incl.cut=1,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE,
                          dcc = FALSE)
ttout3enhanced
# turn remainder rows to 0
ttout3enhanced$tt['3', 'OUT'] <- 0
ttout3enhanced$tt['7', 'OUT'] <- 0
ttout3enhanced$tt['19', 'OUT'] <- 0
ttout3enhanced$tt['23', 'OUT'] <- 0

#enhanced truth table
ttout3enhanced  
#write.csv(ttout3enhanced$tt, "ttout3enhanced.csv")

#enhanced parsimonious solution
psttout3enhanced <- minimize(ttout3enhanced, include="?", details=TRUE, show.cases=TRUE, row.dom=TRUE, all.sol=TRUE, use.tilde=TRUE)
psttout3enhanced 

psttout3 # for comparison

#enhanced intermediate solution
isttout3enhanced <- minimize(ttout3enhanced, include="?", dir.exp="-,-,1,0,1", details=TRUE, show.cases=TRUE, row.dom=TRUE,
                         all.sol=TRUE, use.tilde=TRUE)
isttout3enhanced  ######### FINAL SOLUtION FOR THE OUTCOME "removal from office" #########
isttout3enhanced$i.sol$C1P1$EC 
isttout3enhanced$i.sol$C1P1$DC

isttout3 # original version for comparison

```

# Plot enhanced intermediate solution (Figure 1, Figure 2, Appendix 6 and Appendix 7)
```{r}

## Outcome: out_lower_chamber_vote (see Appendix 6)

# 1 ~elite_ide_polar*mass_mobilization*~leg_shield +
# 2 ~elite_aff_polar*mass_mobilization*corrupt_scandal +
# 3 ~elite_ide_polar*elite_aff_polar*mass_mobilization*~corrupt_scandal +
# 4 ~elite_ide_polar*elite_aff_polar*~leg_shield*corrupt_scandal +
# 5 elite_ide_polar*~elite_aff_polar*~leg_shield*corrupt_scandal +
# 6 elite_ide_polar*elite_aff_polar*~mass_mobilization*~leg_shield*~corrupt_scandal

df_set_out1$term_out2_1 <- as.numeric(fuzzyand(1-df_set_out1$elite_ide_polar, df_set_out1$mass_mobilization, 1-df_set_out1$leg_shield))
df_set_out1$term_out2_2 <- as.numeric(fuzzyand(1-df_set_out1$elite_aff_polar, df_set_out1$mass_mobilization, df_set_out1$corrupt_scandal))
df_set_out1$term_out2_3 <- as.numeric(fuzzyand(1-df_set_out1$elite_ide_polar, df_set_out1$elite_aff_polar, df_set_out1$mass_mobilization, 1-df_set_out1$corrupt_scandal))
df_set_out1$term_out2_4 <- as.numeric(fuzzyand(1-df_set_out1$elite_ide_polar, df_set_out1$elite_aff_polar, 1-df_set_out1$leg_shield, df_set_out1$corrupt_scandal))
df_set_out1$term_out2_5 <- as.numeric(fuzzyand(df_set_out1$elite_ide_polar, 1-df_set_out1$elite_aff_polar, 1-df_set_out1$leg_shield, df_set_out1$corrupt_scandal))
df_set_out1$term_out2_6 <- as.numeric(fuzzyand(df_set_out1$elite_ide_polar, df_set_out1$elite_aff_polar, 1-df_set_out1$mass_mobilization, 1-df_set_out1$leg_shield, 1-df_set_out1$corrupt_scandal))

Plotterm_out2_1 <- XYplot(df_set_out1$term_out2_1, df_set_out1$out_lower_chamber_vote, mguides = TRUE, jitter = TRUE, xlab = "term_out2_ 1", ylab = "outcome", cex = 0.5, clabels = rownames(df_set_out1)) 
Plotterm_out2_2 <- XYplot(df_set_out1$term_out2_2, df_set_out1$out_lower_chamber_vote, mguides = TRUE, jitter = TRUE, xlab = "term_out2_ 2", ylab = "outcome", cex = 0.5, clabels = rownames(df_set_out1))
Plotterm_out2_3 <- XYplot(df_set_out1$term_out2_3, df_set_out1$out_lower_chamber_vote, mguides = TRUE, jitter = TRUE, xlab = "term_out2_ 3", ylab = "outcome", cex = 0.5, clabels = rownames(df_set_out1))
Plotterm_out2_4 <- XYplot(df_set_out1$term_out2_4, df_set_out1$out_lower_chamber_vote, mguides = TRUE, jitter = TRUE, xlab = "term_out2_ 1", ylab = "outcome", cex = 0.5, clabels = rownames(df_set_out1)) 
Plotterm_out2_5 <- XYplot(df_set_out1$term_out2_5, df_set_out1$out_lower_chamber_vote, mguides = TRUE, jitter = TRUE, xlab = "term_out2_ 1", ylab = "outcome", cex = 0.5, clabels = rownames(df_set_out1)) 
Plotterm_out2_6 <- XYplot(df_set_out1$term_out2_6, df_set_out1$out_lower_chamber_vote, mguides = TRUE, jitter = TRUE, xlab = "term_out2_ 1", ylab = "outcome", cex = 0.5, clabels = rownames(df_set_out1)) 

attach(df_set_out1)
df_set_out1$solution_out2 <- fuzzyor(df_set_out1$term_out2_1, df_set_out1$term_out2_2, df_set_out1$term_out2_3, df_set_out1$term_out2_4, df_set_out1$term_out2_5, df_set_out1$term_out2_6)
detach(df_set_out1)

Plotsol_out2 <- XYplot(df_set_out1$solution_out2, df_set_out1$out_lower_chamber_vote, mguides = TRUE, jitter = TRUE, xlab = "Enhanced intermediate solution (Removal from Office)", ylab = "outcome", cex = 0.6, clabels = rownames(df_set_out1))

# stylized full solution plot (optimized for the printed version) (Figure 1)
IS_QCA_out2 <- ggplot(df_set_out1, mapping = aes(x = solution_out2, y = out_lower_chamber_vote)) +
  geom_jitter(width = 0.002, height = 0.002) + #adding jitter to points
  ggrepel::geom_text_repel(aes(label=row.names(df_set_out1)), max.overlaps = Inf, size = 4.5, box.padding = unit(0.8, "lines"), point.padding = unit(0.5, "lines"), segment.color = 'grey50') + #adding labels to points
  geom_segment(aes(x = -0.08, y = -0.11, xend = 1.08, yend = 1.1), linetype = "dotted") +
  geom_hline(yintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  geom_vline (xintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  theme(plot.margin = margin(0.6,0.6,0.6,0.6, "cm")) + ##setting wider plot margins for better viewing
  ggtitle("Sufficiency relation") + #Adding and customising title and axis labels
  xlab("Intermediate Solution")  + ylab("Outcome: Lower Chamber Vote") +
  scale_x_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) + #Customising breaks and axis ticks
  scale_y_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) +
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5, face="bold", size=14, lineheight=1.2), # centering the title and making it bold and larger 
        axis.title.x = element_text(size = 14), # increase x-axis label size
        axis.title.y = element_text(size = 14)) + # increase y-axis label size) 
  annotate("text", x = 0.24, y = 1.08, label = "Inclusion: 0.95  Coverage:  0.655  PRI: 0.95", hjust = 0, size = 5) #adding the information on the right
IS_QCA_out2

ggsave("IS_QCA_out2.png",
       width = 47,
       height = 28,
       units = "cm")

# stylized term 1 (Appendix: Figure 6F)
IS_QCA_out2_term1 <- ggplot(df_set_out1, mapping = aes(x = term_out2_1, y = out_lower_chamber_vote)) +
  geom_jitter(width = 0.002, height = 0.002) + #adding jitter to points
  ggrepel::geom_text_repel(aes(label=row.names(df_set_out1)), max.overlaps = Inf, size = 4.5, box.padding = unit(0.8, "lines"), point.padding = unit(0.5, "lines"), segment.color = 'grey50') + #adding labels to points
  geom_segment(aes(x = -0.08, y = -0.11, xend = 1.08, yend = 1.11), linetype = "dotted") +
  geom_hline(yintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  geom_vline (xintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  theme(plot.margin = margin(0.6,0.6,0.6,0.6, "cm")) + ##setting wider plot margins for better viewing
  ggtitle("Sufficiency relation") + #Adding and customising title and axis labels
  xlab("~elite_ide_polar*mass_mobilization*~leg_shield")  + ylab("Outcome: Lower Chamber Vote") +
  scale_x_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) + #Customising breaks and axis ticks
  scale_y_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) +
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5, face="bold", size=14, lineheight=1.2), # centering the title and making it bold and larger 
        axis.title.x = element_text(size = 14), # increase x-axis label size
        axis.title.y = element_text(size = 14)) + # increase y-axis label size) 
  annotate("text", x = 0.51, y = 1.1, label = "Inclusion: 1.0  Coverage:  0.069  PRI: 1.0", hjust = 0, size = 5) #adding the information on the right
IS_QCA_out2_term1

ggsave("IS_QCA_out2_term1.png",
       width = 47,
       height = 28,
       units = "cm")

# stylized term 2 (Appendix: Figure 6E)
IS_QCA_out2_term2 <- ggplot(df_set_out1, mapping = aes(x = term_out2_2, y = out_lower_chamber_vote)) +
  geom_jitter(width = 0.002, height = 0.002) + #adding jitter to points
  ggrepel::geom_text_repel(aes(label=row.names(df_set_out1)), max.overlaps = Inf, size = 4.5, box.padding = unit(0.8, "lines"), point.padding = unit(0.5, "lines"), segment.color = 'grey50') + #adding labels to points
  geom_segment(aes(x = -0.08, y = -0.11, xend = 1.08, yend = 1.11), linetype = "dotted") +
  geom_hline(yintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  geom_vline (xintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  theme(plot.margin = margin(0.6,0.6,0.6,0.6, "cm")) + ##setting wider plot margins for better viewing
  ggtitle("Sufficiency relation") + #Adding and customising title and axis labels
  xlab("~elite_aff_polar*mass_mobilization*corrupt_scandal")  + ylab("Outcome: Lower Chamber Vote") +
  scale_x_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) + #Customising breaks and axis ticks
  scale_y_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) +
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5, face="bold", size=14, lineheight=1.2), # centering the title and making it bold and larger 
        axis.title.x = element_text(size = 14), # increase x-axis label size
        axis.title.y = element_text(size = 14)) + # increase y-axis label size) 
  annotate("text", x = 0.51, y = 1.1, label = "Inclusion: 1.0  Coverage:  0.138  PRI: 1.0", hjust = 0, size = 5) #adding the information on the right
IS_QCA_out2_term2

ggsave("IS_QCA_out2_term2.png",
       width = 47,
       height = 28,
       units = "cm")

# stylized term 3 (Appendix: Figure 6D)
IS_QCA_out2_term3 <- ggplot(df_set_out1, mapping = aes(x = term_out2_3, y = out_lower_chamber_vote)) +
  geom_jitter(width = 0.002, height = 0.002) + #adding jitter to points
  ggrepel::geom_text_repel(aes(label=row.names(df_set_out1)), max.overlaps = Inf, size = 4.5, box.padding = unit(0.8, "lines"), point.padding = unit(0.5, "lines"), segment.color = 'grey50') + #adding labels to points
  geom_segment(aes(x = -0.08, y = -0.11, xend = 1.08, yend = 1.11), linetype = "dotted") +
  geom_hline(yintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  geom_vline (xintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  theme(plot.margin = margin(0.6,0.6,0.6,0.6, "cm")) + ##setting wider plot margins for better viewing
  ggtitle("Sufficiency relation") + #Adding and customising title and axis labels
  xlab("~elite_ide_polar*elite_aff_polar*mass_mobilization*~corrupt_scandal")  + ylab("Outcome: Lower Chamber Vote") +
  scale_x_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) + #Customising breaks and axis ticks
  scale_y_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) +
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5, face="bold", size=14, lineheight=1.2), # centering the title and making it bold and larger 
        axis.title.x = element_text(size = 14), # increase x-axis label size
        axis.title.y = element_text(size = 14)) + # increase y-axis label size) 
  annotate("text", x = 0.51, y = 1.1, label = "Inclusion: 1.0  Coverage:  0.069  PRI: 1.0", hjust = 0, size = 5) #adding the information on the right
IS_QCA_out2_term3

ggsave("IS_QCA_out2_term3.png",
       width = 47,
       height = 28,
       units = "cm")

# stylized term 4 (Appendix: Figure 6C)
IS_QCA_out2_term4 <- ggplot(df_set_out1, mapping = aes(x = term_out2_4, y = out_lower_chamber_vote)) +
  geom_jitter(width = 0.002, height = 0.002) + #adding jitter to points
  ggrepel::geom_text_repel(aes(label=row.names(df_set_out1)), max.overlaps = Inf, size = 4.5, box.padding = unit(0.8, "lines"), point.padding = unit(0.5, "lines"), segment.color = 'grey50') + #adding labels to points
  geom_segment(aes(x = -0.08, y = -0.11, xend = 1.08, yend = 1.11), linetype = "dotted") +
  geom_hline(yintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  geom_vline (xintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  theme(plot.margin = margin(0.6,0.6,0.6,0.6, "cm")) + ##setting wider plot margins for better viewing
  ggtitle("Sufficiency relation") + #Adding and customising title and axis labels
  xlab("~elite_ide_polar*elite_aff_polar*~leg_shield*corrupt_scandal")  + ylab("Outcome: Lower Chamber Vote") +
  scale_x_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) + #Customising breaks and axis ticks
  scale_y_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) +
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5, face="bold", size=14, lineheight=1.2), # centering the title and making it bold and larger 
        axis.title.x = element_text(size = 14), # increase x-axis label size
        axis.title.y = element_text(size = 14)) + # increase y-axis label size) 
  annotate("text", x = 0.51, y = 1.1, label = "Inclusion: 1.0  Coverage:  0.034  PRI: 1.0", hjust = 0, size = 5) #adding the information on the right
IS_QCA_out2_term4

ggsave("IS_QCA_out2_term4.png",
       width = 47,
       height = 28,
       units = "cm")

# stylized term 5 (Appendix: Figure 6B)
IS_QCA_out2_term5 <- ggplot(df_set_out1, mapping = aes(x = term_out2_5, y = out_lower_chamber_vote)) +
  geom_jitter(width = 0.002, height = 0.002) + #adding jitter to points
  ggrepel::geom_text_repel(aes(label=row.names(df_set_out1)), max.overlaps = Inf, size = 4.5, box.padding = unit(0.8, "lines"), point.padding = unit(0.5, "lines"), segment.color = 'grey50') + #adding labels to points
  geom_segment(aes(x = -0.08, y = -0.11, xend = 1.08, yend = 1.11), linetype = "dotted") +
  geom_hline(yintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  geom_vline (xintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  theme(plot.margin = margin(0.6,0.6,0.6,0.6, "cm")) + ##setting wider plot margins for better viewing
  ggtitle("Sufficiency relation") + #Adding and customising title and axis labels
  xlab("elite_ide_polar*~elite_aff_polar*~leg_shield*corrupt_scandal")  + ylab("Outcome: Lower Chamber Vote") +
  scale_x_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) + #Customising breaks and axis ticks
  scale_y_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) +
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5, face="bold", size=14, lineheight=1.2), # centering the title and making it bold and larger 
        axis.title.x = element_text(size = 14), # increase x-axis label size
        axis.title.y = element_text(size = 14)) + # increase y-axis label size) 
  annotate("text", x = 0.51, y = 1.1, label = "Inclusion: 1.0  Coverage: 0.103  PRI: 1.0", hjust = 0, size = 5) #adding the information on the right
IS_QCA_out2_term5

ggsave("IS_QCA_out2_term5.png",
       width = 47,
       height = 28,
       units = "cm")

# stylized term 6 (Appendix: Figure 6A)
IS_QCA_out2_term6 <- ggplot(df_set_out1, mapping = aes(x = term_out2_6, y = out_lower_chamber_vote)) +
  geom_jitter(width = 0.002, height = 0.002) + #adding jitter to points
  ggrepel::geom_text_repel(aes(label=row.names(df_set_out1)), max.overlaps = Inf, size = 4.5, box.padding = unit(0.8, "lines"), point.padding = unit(0.5, "lines"), segment.color = 'grey50') + #adding labels to points
  geom_segment(aes(x = -0.08, y = -0.11, xend = 1.08, yend = 1.11), linetype = "dotted") +
  geom_hline(yintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  geom_vline (xintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  theme(plot.margin = margin(0.6,0.6,0.6,0.6, "cm")) + ##setting wider plot margins for better viewing
  ggtitle("Sufficiency relation") + #Adding and customising title and axis labels
  xlab("elite_ide_polar*elite_aff_polar*~mass_mobilization*~leg_shield*~corrupt_scandal")  + ylab("Outcome: Lower Chamber Vote") +
  scale_x_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) + #Customising breaks and axis ticks
  scale_y_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) +
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5, face="bold", size=14, lineheight=1.2), # centering the title and making it bold and larger 
        axis.title.x = element_text(size = 14), # increase x-axis label size
        axis.title.y = element_text(size = 14)) + # increase y-axis label size) 
  annotate("text", x = 0.51, y = 1.1, label = "Inclusion: 0.9  Coverage: 0.310  PRI: 0.9", hjust = 0, size = 5) #adding the information on the right
IS_QCA_out2_term6

ggsave("IS_QCA_out2_term6.png",
       width = 47,
       height = 28,
       units = "cm")



## Outcome: out_removal

# 1 ~elite_ide_polar*elite_aff_polar*mass_mobilization +
# 2 elite_aff_polar*~leg_shield*corrupt_scandal +
# 3 ~elite_aff_polar*mass_mobilization*leg_shield*corrupt_scandal

df_set_out2$term_out3_1 <- as.numeric(fuzzyand(1-df_set_out2$elite_ide_polar, df_set_out2$elite_aff_polar, df_set_out2$mass_mobilization))
df_set_out2$term_out3_2 <- as.numeric(fuzzyand(df_set_out2$elite_aff_polar, 1-df_set_out2$leg_shield, df_set_out2$corrupt_scandal))
df_set_out2$term_out3_3 <- as.numeric(fuzzyand(1-df_set_out2$elite_aff_polar, df_set_out2$mass_mobilization, df_set_out2$leg_shield, df_set_out2$corrupt_scandal))

Plotterm_out3_1 <- XYplot(df_set_out2$term_out3_1, df_set_out2$out_removal, mguides = TRUE, jitter = TRUE, xlab = "term_out3_ 1", ylab = "outcome", cex = 0.5, clabels = rownames(df_set_out2)) 
Plotterm_out3_2 <- XYplot(df_set_out2$term_out3_2, df_set_out2$out_removal, mguides = TRUE, jitter = TRUE, xlab = "term_out3_ 2", ylab = "outcome", cex = 0.5, clabels = rownames(df_set_out2))
Plotterm_out3_3 <- XYplot(df_set_out2$term_out3_3, df_set_out2$out_removal, mguides = TRUE, jitter = TRUE, xlab = "term_out3_ 3", ylab = "outcome", cex = 0.5, clabels = rownames(df_set_out2))

attach(df_set_out2)
df_set_out2$solution_out3 <- fuzzyor(df_set_out2$term_out3_1, df_set_out2$term_out3_2, df_set_out2$term_out3_3)
detach(df_set_out2)

Plotsol_out3 <- XYplot(df_set_out2$solution_out3, df_set_out2$out_removal, mguides = TRUE, jitter = TRUE, xlab = "Enhanced intermediate solution (Removal from Office)", ylab = "outcome", cex = 0.6, clabels = rownames(df_set_out2))

# stylized full solution plot (optimized for the printed version) (Figure 2)
IS_QCA_out3 <- ggplot(df_set_out2, mapping = aes(x = solution_out3, y = out_removal)) +
  geom_jitter(width = 0.002, height = 0.002) + #adding jitter to points
  ggrepel::geom_text_repel(aes(label=row.names(df_set_out2)), max.overlaps = Inf, size = 4.5, box.padding = unit(0.8, "lines"), point.padding = unit(0.5, "lines"), segment.color = 'grey50') + #adding labels to points
  geom_segment(aes(x = -0.08, y = -0.11, xend = 1.08, yend = 1.11), linetype = "dotted") +
  geom_hline(yintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  geom_vline (xintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  theme(plot.margin = margin(0.6,0.6,0.6,0.6, "cm")) + ##setting wider plot margins for better viewing
  ggtitle("Sufficiency relation") + #Adding and customising title and axis labels
  xlab("Intermediate Solution")  + ylab("Outcome: Removal from Office") +
  scale_x_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) + #Customising breaks and axis ticks
  scale_y_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) +
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5, face="bold", size=14, lineheight=1.2), # centering the title and making it bold and larger 
        axis.title.x = element_text(size = 14), # increase x-axis label size
        axis.title.y = element_text(size = 14)) + # increase y-axis label size) 
  annotate("text", x = 0.25, y = 1.09, label = "Inclusion: 1.0  Coverage:  0.75  PRI: 1.0", hjust = 0, size = 5) #adding the information on the right
IS_QCA_out3

ggsave("IS_QCA_out3.png",
       width = 47,
       height = 28,
       units = "cm")

# stylized term 1 (Appendix: Figure 7B)
IS_QCA_out3_term1 <- ggplot(df_set_out2, mapping = aes(x = term_out3_1, y = out_removal)) +
  geom_jitter(width = 0.002, height = 0.002) + #adding jitter to points
  ggrepel::geom_text_repel(aes(label=row.names(df_set_out2)), max.overlaps = Inf, size = 4.5, box.padding = unit(0.8, "lines"), point.padding = unit(0.5, "lines"), segment.color = 'grey50') + #adding labels to points
  geom_segment(aes(x = -0.08, y = -0.11, xend = 1.08, yend = 1.11), linetype = "dotted") +
  geom_hline(yintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  geom_vline (xintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  theme(plot.margin = margin(0.6,0.6,0.6,0.6, "cm")) + ##setting wider plot margins for better viewing
  ggtitle("Sufficiency relation") + #Adding and customising title and axis labels
  xlab("~elite_ide_polar*elite_aff_polar*mass_mobilization")  + ylab("Outcome: Removal from Office") +
  scale_x_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) + #Customising breaks and axis ticks
  scale_y_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) +
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5, face="bold", size=14, lineheight=1.2), # centering the title and making it bold and larger 
        axis.title.x = element_text(size = 14), # increase x-axis label size
        axis.title.y = element_text(size = 14)) + # increase y-axis label size) 
  annotate("text", x = 0.51, y = 1.1, label = "Inclusion: 1.0  Coverage:  0.1  PRI: 1.0", hjust = 0, size = 5) #adding the information on the right
IS_QCA_out3_term1

ggsave("IS_QCA_out3_term1.png",
       width = 47,
       height = 28,
       units = "cm")

# stylized term 2 (Appendix: Figure 7A)
IS_QCA_out3_term2 <- ggplot(df_set_out2, mapping = aes(x = term_out3_2, y = out_removal)) +
  geom_jitter(width = 0.002, height = 0.002) + #adding jitter to points
  ggrepel::geom_text_repel(aes(label=row.names(df_set_out2)), max.overlaps = Inf, size = 4.5, box.padding = unit(0.8, "lines"), point.padding = unit(0.5, "lines"), segment.color = 'grey50') + #adding labels to points
  geom_segment(aes(x = -0.08, y = -0.11, xend = 1.08, yend = 1.11), linetype = "dotted") +
  geom_hline(yintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  geom_vline (xintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  theme(plot.margin = margin(0.6,0.6,0.6,0.6, "cm")) + ##setting wider plot margins for better viewing
  ggtitle("Sufficiency relation") + #Adding and customising title and axis labels
  xlab("elite_aff_polar*~leg_shield*corrupt_scandal")  + ylab("Outcome: Removal from Office") +
  scale_x_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) + #Customising breaks and axis ticks
  scale_y_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) +
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5, face="bold", size=14, lineheight=1.2), # centering the title and making it bold and larger 
        axis.title.x = element_text(size = 14), # increase x-axis label size
        axis.title.y = element_text(size = 14)) + # increase y-axis label size) 
  annotate("text", x = 0.51, y = 1.1, label = "Inclusion: 1.0  Coverage:  0.5  PRI: 1.0", hjust = 0, size = 5) #adding the information on the right
IS_QCA_out3_term2

ggsave("IS_QCA_out3_term2.png",
       width = 47,
       height = 28,
       units = "cm")

# stylized term 3 (Appendix: Figure 7C)
IS_QCA_out3_term3 <- ggplot(df_set_out2, mapping = aes(x = term_out3_3, y = out_removal)) +
  geom_jitter(width = 0.002, height = 0.002) + #adding jitter to points
  ggrepel::geom_text_repel(aes(label=row.names(df_set_out2)), max.overlaps = Inf, size = 4.5, box.padding = unit(0.8, "lines"), point.padding = unit(0.5, "lines"), segment.color = 'grey50') + #adding labels to points
  geom_segment(aes(x = -0.08, y = -0.11, xend = 1.08, yend = 1.11), linetype = "dotted") +
  geom_hline(yintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  geom_vline (xintercept = 0.5, linetype = "dashed", alpha = 0.5) +
  theme(plot.margin = margin(0.6,0.6,0.6,0.6, "cm")) + ##setting wider plot margins for better viewing
  ggtitle("Sufficiency relation") + #Adding and customising title and axis labels
  xlab("~elite_aff_polar*mass_mobilization*leg_shield*corrupt_scandal")  + ylab("Outcome: Removal from Office") +
  scale_x_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) + #Customising breaks and axis ticks
  scale_y_continuous(expand = c(0, 0), breaks = seq(0, 1, 0.1)) +
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5, face="bold", size=14, lineheight=1.2), # centering the title and making it bold and larger 
        axis.title.x = element_text(size = 14), # increase x-axis label size
        axis.title.y = element_text(size = 14)) + # increase y-axis label size) 
  annotate("text", x = 0.51, y = 1.1, label = "Inclusion: 1.0  Coverage:  0.15  PRI: 1.0", hjust = 0, size = 5) #adding the information on the right
IS_QCA_out3_term3

ggsave("IS_QCA_out3_term3.png",
       width = 47,
       height = 28,
       units = "cm")

```

# Sufficiency analysis for the non-outcome (standard & enhanced) (Appendix 10)
```{r}

## Outcome: Anti-impeachment lower house vote

ttnoout2 <- truthTable(data=df_set_out1, outcome = "~out_lower_chamber_vote", conditions = "elite_ide_polar, elite_aff_polar, mass_mobilization, leg_shield, corrupt_scandal", 
                          incl.cut=0.9,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE, 
                          dcc = TRUE)
ttnoout2

csttnoout2 <- minimize(ttnoout2, details=TRUE, show.cases=TRUE, row.dom=TRUE, all.sol=TRUE, use.tilde=TRUE)
csttnoout2
csttnoout2$PIchart

psttnoout2 <- minimize(ttnoout2, include="?", details=TRUE, show.cases=TRUE, row.dom=TRUE, all.sol=TRUE, use.tilde=TRUE)
psttnoout2
psttnoout2$PIchart 
psttnoout2$SA

isttnoout2 <- minimize(ttnoout2, include="?", dir.exp="-,-,0,1,0", details=TRUE, show.cases=TRUE, row.dom=TRUE, 
                         all.sol=TRUE, use.tilde=TRUE)
isttnoout2
isttnoout2$i.sol$C1P1$EC 
isttnoout2$i.sol$C1P1$DC

#Enhanced standard analysis for the absence of the Outcome ~out 

#1 Contradictions to necessity:
# no necessary conditions for Oout -> no contradictions possible
#2 Implausible assumptions
# No implausible assumptions
#3 Simultaneous subset relations 
SSR <- intersect(rownames(ttout2$tt) [ttout2$tt$OUT==1],
                                rownames(ttnoout2$tt) [ttnoout2$tt$OUT==1])
SSR
# = 0
# -> no simultaneous subset relations can be found
#4 Incoherent counterfactuals
contcountercrisp2_nout2 <- findRows(obj=ttnoout2, type = 2)
contcountercrisp2_nout2
# Incoherent rows: 4 12 16 20 27
#to be able to compare them, a new truth table is constructed
ttnoout2enhanced <- ttnoout2

# turn remainder rows to 0
ttnoout2enhanced$tt['4', 'OUT'] <- 0
ttnoout2enhanced$tt['12', 'OUT'] <- 0
ttnoout2enhanced$tt['16', 'OUT'] <- 0
ttnoout2enhanced$tt['20', 'OUT'] <- 0
ttnoout2enhanced$tt['27', 'OUT'] <- 0

ttnoout2enhanced
#write.csv(ttnoout2enhanced$tt, "ttnoout2nhanced.csv")

#enhanced parsimonious solution
psttnoout2enhanced <- minimize(ttnoout2enhanced, include="?", details=TRUE, show.cases=TRUE, row.dom=TRUE, all.sol=TRUE, use.tilde=TRUE)
psttnoout2enhanced #enhanced version

#enhanced intermediate solution
isttnoout2enhanced <- minimize(ttnoout2enhanced, include="?", dir.exp="-,-,0,1,0", details=TRUE, show.cases=TRUE, row.dom=TRUE,
                         all.sol=TRUE, use.tilde=TRUE) 
isttnoout2enhanced



## Outcome: Non-removal from office

ttnoout3 <- truthTable(data=df_set_out2, outcome = "~out_removal", conditions = "elite_ide_polar, elite_aff_polar, mass_mobilization, leg_shield, corrupt_scandal", 
                          incl.cut=1,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE,
                          dcc = TRUE)
ttnoout3

csttnoout3 <- minimize(ttnoout3, details=TRUE, show.cases=TRUE, row.dom=TRUE, all.sol=TRUE, use.tilde=TRUE)
csttnoout3
csttnoout3$PIchart

psttnoout3 <- minimize(ttnoout3, include="?", details=TRUE, show.cases=TRUE, row.dom=TRUE, all.sol=TRUE, use.tilde=TRUE)
psttnoout3
psttnoout3$PIchart 
psttnoout3$SA

isttnoout3 <- minimize(ttnoout3, include="?", dir.exp="-,-,0,1,0", details=TRUE, show.cases=TRUE, row.dom=TRUE, 
                         all.sol=TRUE, use.tilde=TRUE)
isttnoout3
isttnoout3$i.sol$C1P1$EC 
isttnoout3$i.sol$C1P1$DC

#Enhanced standard analysis for the absence of the Outcome ~out 

#1 Contradictions to necessity:
# no necessary conditions for Out -> no contradictions possible
#2 Implausible assumptions
# No implausible assumptions
#3 Simultaneous subset relations 
SSR <- intersect(rownames(ttout3$tt) [ttout3$tt$OUT==1],
                                rownames(ttnoout3$tt) [ttnoout3$tt$OUT==1])
SSR
# = 0
# -> no simultaneous subset relations can be found
#4 Incoherent counterfactuals
contcountercrisp2_nout3 <- findRows(obj=ttnoout3, type = 2)
contcountercrisp2_nout3
# Incoherent rows: 3  7 19 23
#to be able to compare them, a new truth table is constructed
ttnoout3enhanced <- ttnoout3

# turn remainder rows to 0
ttnoout3enhanced$tt['3', 'OUT'] <- 0
ttnoout3enhanced$tt['7', 'OUT'] <- 0
ttnoout3enhanced$tt['19', 'OUT'] <- 0
ttnoout3enhanced$tt['23', 'OUT'] <- 0
ttnoout3enhanced
#write.csv(ttnoout3enhanced$tt, "ttnoout3enhanced.csv")

#enhanced parsimonious solution
psttnoout3enhanced <- minimize(ttnoout3enhanced, include="?", details=TRUE, show.cases=TRUE, row.dom=TRUE, all.sol=TRUE, use.tilde=TRUE)
psttnoout3enhanced #enhanced version

#enhanced intermediate solution
isttnoout3enhanced <- minimize(ttnoout3enhanced, include="?", dir.exp="-,-,0,1,0", details=TRUE, show.cases=TRUE, row.dom=TRUE,
                         all.sol=TRUE, use.tilde=TRUE) 
isttnoout3enhanced

```

# Robustness tests (Appendix 11)
```{r}
## Raw consistency robustness
# No need for out 3 since consistency = 1.0

rob.inclrange(
  data = df_set_out1,
  step = 0.01,
  max.runs = 10,
  outcome = "out_lower_chamber_vote",
  conditions = c("elite_ide_polar", "elite_aff_polar", "mass_mobilization", "leg_shield", "corrupt_scandal"),
  incl.cut = 0.9,
  n.cut = 1,
  include = "?", 
  dir.exp = c("-","-",1,0,1))
# Raw Consistency T.:  Lower bound  NA Threshold  0.9 Upper bound  0.9 


## Frequency cut-off robustness
rob.ncutrange(
  data = df_set_out2,
  step = 1,
  max.runs = 10,
  outcome = "out_removal",
  conditions = c("elite_ide_polar", "elite_aff_polar", "mass_mobilization", "leg_shield", "corrupt_scandal"),
  incl.cut = 1.0,
  n.cut = 1,
  include = "?", 
  dir.exp = c("-","-",1,0,1))
# N.Cut:  Lower bound  1 Threshold  1 Upper bound  1 

rob.ncutrange(
  data = df_set_out1,
  step = 1,
  max.runs = 10,
  outcome = "out_lower_chamber_vote",
  conditions = c("elite_ide_polar", "elite_aff_polar", "mass_mobilization", "leg_shield", "corrupt_scandal"),
  incl.cut = 0.9,
  n.cut = 1,
  include = "?", 
  dir.exp = c("-","-",1,0,1))
# N.Cut:  Lower bound  1 Threshold  1 Upper bound  1 


## sensitivity range for calibration

# subset raw data
data_raw_out1 <- data_r
data_raw_out2 <- subset(data_r, out_lower_chamber_vote == 1)

# assign president-years as rownames
data_raw_out1 <- data.frame(data_raw_out1[,-1], row.names=data_raw_out1[,1])
data_raw_out2 <- data.frame(data_raw_out2[,-1], row.names=data_raw_out2[,1])


# Original cliberation thresholds
# elite_ide_polar  5
# elite_aff_polar 1.001
# mass_mobilization 10000
# leg_shield 0
# corrupt_scandal 0.5 (binary) # no need for sensitivity range

# Sensitivity range calibration: elite ideological polarization
rob.calibrange(raw.data = data_raw_out2,
               calib.data = df_set_out2, 
               test.cond.raw = "elite_ide_polar",
               test.cond.calib = "elite_ide_polar",
               test.thresholds = c(5),
               type = "crisp",
               step = 0.1,
               max.runs = 40,
               outcome = "out_removal",
               conditions = c("elite_ide_polar", "elite_aff_polar", "mass_mobilization", "leg_shield", "corrupt_scandal"),
               incl.cut = 1.0,
               n.cut = 1,
               include = "?")
# Crossover:  Crossover:  Lower bound  4.9 Threshold  5 Upper bound  5.2 

rob.calibrange(raw.data = data_raw_out1,
               calib.data = df_set_out1, 
               test.cond.raw = "elite_ide_polar",
               test.cond.calib = "elite_ide_polar",
               test.thresholds = c(5), 
               type = "crisp",
               step = 0.1,
               max.runs = 40,
               outcome = "out_lower_chamber_vote",
               conditions = c("elite_ide_polar", "elite_aff_polar", "mass_mobilization", "leg_shield", "corrupt_scandal"),
               incl.cut = 0.9,
               n.cut = 1,
               include = "?")
# Crossover:  Lower bound  4.8 Threshold  5 Upper bound  5.2 

# Sensitivity range calibration: elite affective polarization
data_raw_out2$elite_aff_polar_rev <- 4-data_raw_out2$elite_aff_polar # reversed raw
data_raw_out1$elite_aff_polar_rev <- 4-data_raw_out1$elite_aff_polar # reversed raw

rob.calibrange(raw.data = data_raw_out2,
               calib.data = df_set_out2, 
               test.cond.raw = "elite_aff_polar_rev",
               test.cond.calib = "elite_aff_polar",
               test.thresholds = c(1.001),
               type = "crisp",
               step = 0.1,
               max.runs = 40,
               outcome = "out_removal",
               conditions = c("elite_ide_polar", "elite_aff_polar", "mass_mobilization", "leg_shield", "corrupt_scandal"),
               incl.cut = 1.0,
               n.cut = 1,
               include = "?")
# Crossover:  Lower bound  1.001 Threshold  1.001 Upper bound  1.401 

rob.calibrange(raw.data = data_raw_out1,
               calib.data = df_set_out1, 
               test.cond.raw = "elite_aff_polar_rev",
               test.cond.calib = "elite_aff_polar",
               test.thresholds = c(1.001),
               type = "crisp",
               step = 0.1,
               max.runs = 40,
               outcome = "out_removal",
               conditions = c("elite_ide_polar", "elite_aff_polar", "mass_mobilization", "leg_shield", "corrupt_scandal"),
               incl.cut = 0.9,
               n.cut = 1,
               include = "?")
# Crossover:  Lower bound  1.001 Threshold  1.001 Upper bound  1.401 

# Sensitivity range calibration: mass mooblization
rob.calibrange(raw.data = data_raw_out2,
               calib.data = df_set_out2, 
               test.cond.raw = "mass_mobilization",
               test.cond.calib = "mass_mobilization",
               test.thresholds = c(10000),
               type = "crisp",
               step = 1000,
               max.runs = 40,
               outcome = "out_removal",
               conditions = c("elite_ide_polar", "elite_aff_polar", "mass_mobilization", "leg_shield", "corrupt_scandal"),
               incl.cut = 1.0,
               n.cut = 1,
               include = "?")
# Crossover:  Lower bound  1000 Threshold  10000 Upper bound  NA  

rob.calibrange(raw.data = data_raw_out1,
               calib.data = df_set_out1, 
               test.cond.raw = "mass_mobilization",
               test.cond.calib = "mass_mobilization",
               test.thresholds = c(10000), 
               type = "crisp",
               step = 1000,
               max.runs = 40,
               outcome = "out_lower_chamber_vote",
               conditions = c("elite_ide_polar", "elite_aff_polar", "mass_mobilization", "leg_shield", "corrupt_scandal"),
               incl.cut = 0.9,
               n.cut = 1,
               include = "?")
# Crossover:  Lower bound  6000 Threshold  10000 Upper bound  NA 

# Sensitivity range calibration: legislative shield
## Note: Not working. Error message Threshold value(s) outside the range of x
rob.calibrange(raw.data = data_raw_out2,
               calib.data = df_set_out2, 
               test.cond.raw = "leg_shield",
               test.cond.calib = "leg_shield",
               test.thresholds = c(0.1),
               type = "crisp",
               step = 1,
               max.runs = 40,
               outcome = "out_removal",
               conditions = c("elite_ide_polar", "elite_aff_polar", "mass_mobilization", "leg_shield", "corrupt_scandal"),
               incl.cut = 1.0,
               n.cut = 1,
               include = "?")

rob.calibrange(raw.data = data_raw_out1,
               calib.data = df_set_out1, 
               test.cond.raw = "leg_shield",
               test.cond.calib = "leg_shield",
               test.thresholds = c(0.1), 
               type = "crisp",
               step = 1,
               max.runs = 40,
               outcome = "out_lower_chamber_vote",
               conditions = c("elite_ide_polar", "elite_aff_polar", "mass_mobilization", "leg_shield", "corrupt_scandal"),
               incl.cut = 0.9,
               n.cut = 1,
               include = "?")

#### Fit-oriented and case-oriented robustness

# New original data frame with uppercase for rob.fit to work

df_orig <- data.frame(
  code = data_r$code,
  out_lower_chamber_vote = data_r$out_lower_chamber_vote,
  out_removal = data_r$out_removal,
  elite_ide_polar = calibrate(data_r$elite_ide_polar, 
                              thresholds = 5, 
                              type = "crisp"),
  elite_aff_polar = calibrate(4-data_r$elite_aff_polar, 
                              thresholds = 1.001, 
                              type = "crisp"),
  mass_mobilization = calibrate(data_r$mass_mobilization, 
                                thresholds = 10000, 
                                type = "crisp"),
  leg_shield = calibrate(data_r$leg_shield, 
                         thresholds = 0.1, 
                         type = "crisp"),
  corrupt_scandal = data_r$corrupt_scandal
)

names(df_orig) <- toupper(names(df_orig)) # for the rob.fit below to work

# subset data
df_set_out1_orig <- df_orig
df_set_out2_orig <- subset(df_orig, OUT_LOWER_CHAMBER_VOTE == 1)

# assign president-years as rownames
df_set_out1_orig <- data.frame(df_set_out1_orig[,-1], row.names=df_set_out1_orig[,1])
df_set_out2_orig <- data.frame(df_set_out2_orig[,-1], row.names=df_set_out2_orig[,1])


# Creating original enhanced intermediate solution

## Outcome: Lower house vote

ttout2enhanced_orig <- truthTable(data=df_set_out1_orig, outcome = "OUT_LOWER_CHAMBER_VOTE", conditions = "ELITE_IDE_POLAR, ELITE_AFF_POLAR, MASS_MOBILIZATION, LEG_SHIELD, CORRUPT_SCANDAL", 
                          incl.cut=0.9,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE,
                          dcc = TRUE)
ttout2enhanced_orig

contcountercrisp_out2 <- findRows(obj=ttout2enhanced_orig, type = 2)
contcountercrisp_out2
# Incoherent rows: 4 12 16 20 27

# turn remainder rows to 0
ttout2enhanced_orig$tt['4', 'OUT'] <- 0
ttout2enhanced_orig$tt['12', 'OUT'] <- 0
ttout2enhanced_orig$tt['16', 'OUT'] <- 0
ttout2enhanced_orig$tt['20', 'OUT'] <- 0
ttout2enhanced_orig$tt['27', 'OUT'] <- 0

#enhanced truth table
ttout2enhanced_orig

#enhanced intermediate solution
isttout2enhanced_orig <- minimize(ttout2enhanced_orig, include="?", dir.exp="-,-,1,0,1", details=TRUE, show.cases=TRUE, row.dom=TRUE,
                         all.sol=TRUE, use.tilde=TRUE)
isttout2enhanced_orig 

# Creating the test solution:

# altering consistency
# not necessary given high consistency for both out3 (1.0) and out2 (0.9)

# altering calibration for elite ideological and affective polarization (df_set_out1_rob was created above)
ttout2enhanced_rob <- truthTable(data=df_set_out1_rob, outcome = "OUT_LOWER_CHAMBER_VOTE", conditions = "ELITE_IDE_POLAR, ELITE_AFF_POLAR, MASS_MOBILIZATION, LEG_SHIELD, CORRUPT_SCANDAL", 
                          incl.cut=0.9,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE,
                          dcc = TRUE)
ttout2enhanced_rob

# enhanced robustness intermediate solution
ttnoout3enhanced_rob <- truthTable(data=df_set_out1_rob, outcome = "~OUT_LOWER_CHAMBER_VOTE", conditions = "ELITE_IDE_POLAR, ELITE_AFF_POLAR, MASS_MOBILIZATION, LEG_SHIELD, CORRUPT_SCANDAL", 
                          incl.cut=0.9,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE,
                          dcc = TRUE)
ttnoout3enhanced_rob

# simultaneous subset relations
SSR <- intersect(rownames(ttout2enhanced_rob$tt) [ttout2enhanced_rob$tt$OUT==1],
                                rownames(ttnoout3enhanced_rob$tt) [ttnoout3enhanced_rob$tt$OUT==1])
SSR
# = 0
# -> no simultaneous subset relations can be found

# incoherent counteractuals
contcountercrisp_out3 <- findRows(obj=ttout2enhanced_rob, type = 2)
contcountercrisp_out3
# Incoherent rows: 4  5 12 16 20 27 31

# turn remainder rows to 0
ttout2enhanced_rob$tt['4', 'OUT'] <- 0
ttout2enhanced_rob$tt['5', 'OUT'] <- 0
ttout2enhanced_rob$tt['12', 'OUT'] <- 0
ttout2enhanced_rob$tt['16', 'OUT'] <- 0
ttout2enhanced_rob$tt['20', 'OUT'] <- 0
ttout2enhanced_rob$tt['27', 'OUT'] <- 0
ttout2enhanced_rob$tt['32', 'OUT'] <- 0

#enhanced truth table
ttout2enhanced_rob

#enhanced intermediate solution
isttout2enhanced_rob <- minimize(ttout2enhanced_rob, include="?", dir.exp="-,-,1,0,1", details=TRUE, show.cases=TRUE, row.dom=TRUE,
                         all.sol=TRUE, use.tilde=TRUE)
isttout2enhanced_rob 

# Calculate parameters for the robust core:
rob.corefit(test_sol = isttout2enhanced_rob, initial_sol = isttout2enhanced_orig, outcome = "OUT_LOWER_CHAMBER_VOTE")

# Calculate robustness parameters:
rob_par_out2 <- rob.fit(test_sol = isttout2enhanced_rob, 
              initial_sol = isttout2enhanced_orig,  
              outcome = "OUT_LOWER_CHAMBER_VOTE")
rob_par_out2

# Plotting the initial solution against the test set (Appendix: Figure 11A):
rob.xyplot_out2 <- rob.xyplot(test_sol = isttout2enhanced_rob, 
           initial_sol = isttout2enhanced_orig, 
           outcome = "OUT_LOWER_CHAMBER_VOTE", 
           fontsize = 3.5, 
           jitter=TRUE,
           all_labels = FALSE,
           area_lab=TRUE)
rob.xyplot_out2

ggsave("rob.xyplot_out2.png",
       width = 20,
       height = 16,
       units = "cm")

# Obtaining names of case types and robustness case parameters:
# Note: Beware that the code below takes a while to run
rob.cases(test_sol = isttout2enhanced_rob, 
          initial_sol = isttout2enhanced_orig, 
          outcome = "OUT_LOWER_CHAMBER_VOTE")

# Outcome: Removal from Office

ttout3enhanced_orig <- truthTable(data=df_set_out2_orig, outcome = "OUT_REMOVAL", conditions = "ELITE_IDE_POLAR, ELITE_AFF_POLAR, MASS_MOBILIZATION, LEG_SHIELD, CORRUPT_SCANDAL", 
                          incl.cut=1,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE,
                          dcc = TRUE)

ttout3enhanced_orig

contcountercrisp_out3 <- findRows(obj=ttout3enhanced_orig, type = 2)
contcountercrisp_out3
# Incoherent rows: 3  7 19 23

# turn remainder rows to 0
ttout3enhanced_orig$tt['3', 'OUT'] <- 0
ttout3enhanced_orig$tt['7', 'OUT'] <- 0
ttout3enhanced_orig$tt['19', 'OUT'] <- 0
ttout3enhanced_orig$tt['23', 'OUT'] <- 0

#enhanced truth table
ttout3enhanced_orig

#enhanced intermediate solution
isttout3enhanced_orig <- minimize(ttout3enhanced_orig, include="?", dir.exp="-,-,1,0,1", details=TRUE, show.cases=TRUE, row.dom=TRUE,
                         all.sol=TRUE, use.tilde=TRUE)
isttout3enhanced_orig 

# Creating the test solution:

# altering consistency
# not necessary given high consistency for both out3 (1.0) and out2 (0.9)

# altering calibration for elite ideological and affective polarization
df_rob <- data.frame(
  code = data_r$code,
  out_lower_chamber_vote = data_r$out_lower_chamber_vote,
  out_removal = data_r$out_removal,
  elite_ide_polar = calibrate(data_r$elite_ide_polar, 
                              thresholds = 4.001, # new cross over point
                              type = "crisp"),
  elite_aff_polar = calibrate(4-data_r$elite_aff_polar, # "4-" because it is a reverse scale 
                              thresholds = 1.501, # new cross over point (2.499 in the non-reversed raw data)
                              type = "crisp"),
  mass_mobilization = calibrate(data_r$mass_mobilization, 
                                thresholds = 10000, 
                                type = "crisp"),
  leg_shield = calibrate(data_r$leg_shield, 
                         thresholds = 0.1, 
                         type = "crisp"),
  corrupt_scandal = data_r$corrupt_scandal
)

names(df_rob) <- toupper(names(df_rob)) # for the rob.fit below to work

# subset data
df_set_out1_rob <- df_rob
df_set_out2_rob <- subset(df_rob, OUT_LOWER_CHAMBER_VOTE == 1)

# assign president-years as rownames
df_set_out1_rob <- data.frame(df_set_out1_rob[,-1], row.names=df_set_out1_rob[,1])
df_set_out2_rob <- data.frame(df_set_out2_rob[,-1], row.names=df_set_out2_rob[,1])

ttout3enhanced_rob <- truthTable(data=df_set_out2_rob, outcome = "OUT_REMOVAL", conditions = "ELITE_IDE_POLAR, ELITE_AFF_POLAR, MASS_MOBILIZATION, LEG_SHIELD, CORRUPT_SCANDAL", 
                          incl.cut=1,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE,
                          dcc = TRUE)
ttout3enhanced_rob

# enhanced robustness intermediate solution
ttnoout3enhanced_rob <- truthTable(data=df_set_out2_rob, outcome = "~OUT_REMOVAL", conditions = "ELITE_IDE_POLAR, ELITE_AFF_POLAR, MASS_MOBILIZATION, LEG_SHIELD, CORRUPT_SCANDAL", 
                          incl.cut=1,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE,
                          dcc = TRUE)
ttnoout3enhanced_rob

# simultaneous subset relations
SSR <- intersect(rownames(ttout3enhanced_rob$tt) [ttout3enhanced_rob$tt$OUT==1],
                                rownames(ttnoout3enhanced_rob$tt) [ttnoout3enhanced_rob$tt$OUT==1])
SSR
# = 0
# -> no simultaneous subset relations can be found

# incoherent counteractuals
contcountercrisp_out3 <- findRows(obj=ttout3enhanced_rob, type = 2)
contcountercrisp_out3
# Incoherent rows: 3 19

ttnoout3enhanced_rob$tt['3', 'OUT'] <- 0
ttnoout3enhanced_rob$tt['19', 'OUT'] <- 0

#enhanced intermediate solution
isttout3enhanced_rob <- minimize(ttout3enhanced_rob, include="?", dir.exp="-,-,1,0,1", details=TRUE, show.cases=TRUE, row.dom=TRUE,
                         all.sol=TRUE, use.tilde=TRUE)
isttout3enhanced_rob 

# Calculate parameters for the robust core:
rob.corefit(test_sol = isttout3enhanced_rob, initial_sol = isttout3enhanced_orig, outcome = "OUT_REMOVAL")

# Calculate robustness parameters:
rob_par_out3 <- rob.fit(test_sol = isttout3enhanced_rob, 
              initial_sol = isttout3enhanced_orig,  
              outcome = "OUT_REMOvAL")
rob_par_out3

# Plotting the initial solution against the test set (Appendix: Figure 11B):
rob.xyplot_out3 <- rob.xyplot(test_sol = isttout3enhanced_rob, 
           initial_sol = isttout3enhanced_orig, 
           outcome = "OUT_REMOVAL", 
           fontsize = 3.5, 
           jitter=TRUE,
           all_labels = FALSE,
           area_lab=TRUE)
rob.xyplot_out3

ggsave("rob.xyplot_out3.png",
       width = 20,
       height = 16,
       units = "cm")

# Obtaining names of case types and robustness case parameters:
# Note: Beware that the code below takes a while to run
rob.cases(test_sol = isttout3enhanced_rob, 
          initial_sol = isttout3enhanced_orig, 
          outcome = "OUT_REMOVAL")


```

# Cluster analysis (Figure 3 and Figure 4)
```{r}
# region codes (V-Dem) for e_regionpol_7C
# 1: Eastern Europe (including German Democratic Republic, excluding the Caucasus)
# 2: Latin America and the Caribbean
# 3: The Middle East and North Africa (including Israel and Türkiye, excluding Cyprus)
# 4: Sub-Saharan Africa
# 5: Western Europe and North America (including Cyprus, but excluding German Democratic Republic)
# 6: East Asia and the Pacific
# 7: South and Central Asia (including the Caucasus)


# Outcome: Lower house vote
df_set_out1$code <- rownames(df_set_out1)

CSout2 <- cluster(df_set_out1,
               "~elite_ide_polar*mass_mobilization*~leg_shield + ~elite_aff_polar*mass_mobilization*corrupt_scandal + ~elite_ide_polar*elite_aff_polar*mass_mobilization*~corrupt_scandal + ~elite_ide_polar*elite_aff_polar*~leg_shield*corrupt_scandal + elite_ide_polar*~elite_aff_polar*~leg_shield*corrupt_scandal + elite_ide_polar*elite_aff_polar*~mass_mobilization*~leg_shield*~corrupt_scandal", 
              outcome = "out_lower_chamber_vote",
              unit_id = "code",
              cluster_id = "e_regionpol_7C",
              wicons = FALSE)
CSout2

CSout2_plot <- cluster.plot(cluster.res = CSout2,
             labs = FALSE, # labs TRUE does not work somehow here so I'll create the plot with ggplot
             size = 8,
             angle = 6,
             wicons = FALSE)

## create plot
CSout2style <- data.frame(Cases = c("Africa (5)", "Asia (6)", "Eastern Europe (12)", "Latin America (17)", "North America (3)"), # South and Central Asia excluded because only one case
                      Consistency = c(1, 1, 1, 0.833, 1))

CSout2plotstyle <- ggplot(CSout2style, aes(x = Cases, y = Consistency)) +
  geom_point() +
  geom_hline(yintercept = 0.95, 
             linetype = "dashed", color = "red") +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Clusters", y = "Consistency") +
  theme_bw()

CSout2plotstyle

ggsave("IS_CSout2plotstyle.png",
       width = 20,
       height = 16,
       units = "cm")

# Outcome: Removal from office

# change row_names to a variable in df_set_out2 for the cluster analysis
df_set_out2$code <- rownames(df_set_out2)

CSout3 <- cluster(df_set_out2, "~elite_ide_polar*elite_aff_polar*mass_mobilization + elite_aff_polar*~leg_shield*corrupt_scandal + ~elite_aff_polar*mass_mobilization*leg_shield*corrupt_scandal", 
              outcome = "out_removal",
              unit_id = "code",
              cluster_id = "e_regionpol_7C",
              wicons = FALSE)
CSout3

CSout3_plot <- cluster.plot(cluster.res = CSout3,
             labs = FALSE, # labs TRUE does not work somehow here so I'll create the plot with ggplot
             size = 8,
             angle = 6,
             wicons = FALSE)

## create plot (Figure 4)
CSout3style <- data.frame(Cases = c("Africa (2)", "Asia (4)", "Eastern Europe (8)", "Latin America (11)", "North America (3)"), # South and Central Asia excluded because only one case
                      Consistency = c(1, 1, 1, 1, 1))

CSout3plotstyle <- ggplot(CSout3style, aes(x = Cases, y = Consistency)) +
  geom_point() +
  geom_hline(yintercept = 1, 
             linetype = "dashed", color = "red") +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Clusters", y = "Consistency") +
  theme_bw()

CSout3plotstyle

ggsave("IS_CSout3plotstyle.png",
       width = 20,
       height = 16,
       units = "cm")


```

# Theory evaluation (Appendix 12)
```{r}

## Theory evaluation for out_lower_chamber_vote

theory_out2 <- "ELITE_IDE_POLAR*CORRUPT_SCANDAL + ELITE_IDE_POLAR*CORRUPT_SCANDAL*MASS_MOBILIZATION + ELITE_AFF_POLAR + ELITE_AFF_POLAR*~CORRUPT_SCANDAL + MASS_MOBILIZATION + CORRUPT_SCANDAL*~LEG_SHIELD"

# Perform theory evaluation with fit:
# Note: It takes sometime to run
TEV_out2 <- theory.evaluation(theory = theory_out2, empirics = isttout2enhanced_orig, outcome = "OUT_LOWER_CHAMBER_VOTE", print.fit = TRUE, print.data = TRUE)

TEV_out2

## Theory evaluation for out_removal

# Summarize the theory in Boolean terms:
theory_out3 <- "ELITE_IDE_POLAR*CORRUPT_SCANDAL + ELITE_IDE_POLAR*CORRUPT_SCANDAL*MASS_MOBILIZATION + ELITE_AFF_POLAR + ELITE_AFF_POLAR*~CORRUPT_SCANDAL + MASS_MOBILIZATION + CORRUPT_SCANDAL*~LEG_SHIELD"

# Perform theory evaluation with fit:
# Note: It takes sometime to run
TEV_out3 <- theory.evaluation(theory = theory_out3, empirics = isttout3enhanced_orig,
                         outcome = "OUT_REMOVAL", print.fit = TRUE, print.data = TRUE)

TEV_out3

```

# Cross-varlidation of ideological polarization data
```{r}
# DW-Nominate and V-Party

# load packages
library(tidyverse)
library(dplyr)
library(data.table)

# load v-party dataset

library(vdemdata)
vdemdata::vparty -> vparty
head(vparty)

# load dw-nominate congress

dw_nom_cong <- read.csv(file = "dw_nominate_usa.csv", header=TRUE)
head(dw_nom_cong)

dw_clean <- dw_nom_cong %>%
  filter(congress >= 91, congress <= 115) %>%
  filter(
    party_name %in% c("Democrat", "Republican"),
    chamber != "President") %>% 
  mutate(
    year = 1970 + 2 * (congress - 91),
    year_end = year + 1
  ) %>% 
  mutate(
    party_name = case_when(
      party_name == "Democrat"   ~ "Democratic Party",
      party_name == "Republican" ~ "Republican Party",
      TRUE ~ party_name
    ))
head(dw_clean)

vparty_usa <- vparty %>% 
  select(
    country_text_id, year, v2pariglef,
    party_name = v2paenname
  ) %>% 
  filter(
    country_text_id == "USA",
    year >= 1970
  )
head(vparty_usa)

# Join datasets
dw_vparty <- dw_clean %>% left_join(vparty_usa, by = c("year", "party_name"))
colnames(dw_vparty)

# Correlations by party (raw data)
cor(dw_vparty$nominate_dim1_mean, dw_vparty$v2pariglef, method = "pearson")
# 0.981792
cor(dw_vparty$nominate_dim1_mean, dw_vparty$v2pariglef, method = "spearman")
# 0.8426763

# Correlations by party (rescaled data)
rescale01 <- function(x) (x - min(x, na.rm = TRUE)) /
                          (max(x, na.rm = TRUE) - min(x, na.rm = TRUE))

dw_vparty <- dw_vparty %>%
  mutate(
    nominate_dim1_mean_std = rescale01(nominate_dim1_mean),
    v2pariglef_std = rescale01(v2pariglef)
  )

cor(dw_vparty$nominate_dim1_mean_std, dw_vparty$v2pariglef_std, method = "pearson")
# 0.981792
cor(dw_vparty$nominate_dim1_mean_std, dw_vparty$v2pariglef_std, method = "spearman")
# 0.8426763

# Manifesto Project and V-Party

# Load manifesto project data

manifesto_raw <- read.csv(file = "MPDataset.csv", header=TRUE)
colnames(manifesto_raw)

manifesto_clean <- manifesto_raw %>% 
  select(country_name = countryname, date, rile, planeco, markeco, welfare, party_name = partyname, party_id = party) %>% 
  mutate(year = date %/% 100) %>% 
  select(-date)

# load vparty dataset
Vparty_clean <- vparty %>% 
  select(country_name, year, party_name = v2paenname,  v2pariglef)

# join datasets
manifesto_vparty <- manifesto_clean %>% left_join(Vparty_clean, by = c("country_name", "year", "party_name")) %>%  drop_na() 

# CHES and V-Party

# load CHES data
ches_raw <- read.csv(file = "CHES_dataset_means.csv", header=TRUE)

ches_clean <- ches_raw %>% select(year, party_id = cmp_id, lrgen, lrecon)

ches_manifesto <- manifesto_clean %>% left_join(ches_clean, by = c("party_id", "year")) %>%  drop_na() 

# join datasets
ches_vparty <- ches_manifesto %>% left_join(Vparty_clean, by = c("country_name", "year", "party_name")) %>%  drop_na() 

# Correlations 
cor(ches_vparty$lrecon, ches_vparty$v2pariglef, method = "pearson")
# 0.8230311
cor(ches_vparty$lrecon, ches_vparty$v2pariglef, method = "spearman")
# 0.8268535

```

# Robustness analysis using simulated data (Appendix 13)
```{r}
# load data
data_impeach_sim <- read.csv("set_data_impeach_sim.csv")

attach(data_impeach_sim)  #attach is used to save a bit time, you do not have to write data_raw$year and so on. 

df_set_out2_sim <- data.frame(code, 
                          out_removal_sim = out_removal, 
                          mass_mobilization_sim = mass_mobilization, 
                          leg_shield_sim = leg_shield, 
                          elite_ide_polar_sim = elite_ide_polar, 
                          elite_aff_polar_sim = elite_aff_polar, 
                          corrupt_scandal_sim = corrupt_scandal)

detach(data_impeach_sim)

# assign president-years as rownames
df_set_out2_sim <- data.frame(df_set_out2_sim[,-1], row.names=df_set_out2_sim[,1])
  
# check skewness
check_skewness(df_set_out2_sim)

# --- Necessity analysis --- #

# Out3

cons_out3_sim <- subset(df_set_out2_sim, select = c("elite_ide_polar_sim", "elite_aff_polar_sim", "mass_mobilization_sim", "leg_shield_sim", "corrupt_scandal_sim"))

pof(cons_out3_sim, "out_removal_sim", df_set_out2_sim, relation = "nec")
pof(1-cons_out3_sim, "out_removal_sim", df_set_out2_sim, relation = "nec")
pof(cons_out3_sim, "~out_removal_sim", df_set_out2_sim, relation = "nec")
pof(1-cons_out3_sim, "~out_removal_sim", df_set_out2_sim, relation = "nec")
#no single necessary condition could be found, neither for the outcome, nor for the non-outcome

#check for SUIN condition
SUIN_out3_sim <- superSubset(df_set_out2_sim, outcome = "out_removal_sim", conditions = "elite_ide_polar_sim, elite_aff_polar_sim, mass_mobilization_sim, leg_shield_sim, corrupt_scandal_sim", relation = "nec", incl.cut = 0.9, ron.cut = 0.5, cov.cut = 0.6)
SUIN_out3_sim

#absence of outcome
SUIN_noout3_sim <- superSubset(df_set_out2_sim, outcome = "~out_removal_sim", conditions = "elite_ide_polar_sim, elite_aff_polar_sim, mass_mobilization_sim, leg_shield_sim, corrupt_scandal_sim", relation = "nec", incl.cut = 0.9, ron.cut = 0.5, cov.cut = 0.6)
SUIN_noout3_sim

# --- Sufficiency analysis --- #

# truth tables for out_removal_sim 
ttout3_sim <- truthTable(data=df_set_out2_sim, outcome = "out_removal_sim", conditions = "elite_ide_polar_sim, elite_aff_polar_sim, mass_mobilization_sim, leg_shield_sim, corrupt_scandal_sim", 
                          incl.cut=1,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE,
                          dcc = FALSE)
ttout3_sim

# --- Sufficiency analysis --- #

## out3: Removal

# intermediate solution isttout3_sim (with directional expectations)
isttout3_sim <- minimize(ttout3_sim, include="?", dir.exp="-,-,1,0,1", details=TRUE, show.cases=TRUE, row.dom=TRUE,
                         all.sol=TRUE, use.tilde=TRUE)
isttout3_sim
isttout3_sim$PIchart
isttout3_sim$i.sol$C1P1$EC 
isttout3_sim$i.sol$C1P1$DC

# --- Enhanced intermediate solution --- #

# Steps:
#1 Contradictions to necessity:
# no necessary conditions for the outcome -> no contradictions possible
#2 Implausible assumptions
# There is no possible implausible assumption in my opinion. 
#3 Simultaneous subset relations 
# compare if any of the rows that you included into the solution is also part of a possible solution for the non-outcome
ttnoout3_sim <- truthTable(data=df_set_out2_sim, outcome = "~out_removal_sim", conditions = "elite_ide_polar_sim, elite_aff_polar_sim, mass_mobilization_sim, leg_shield_sim, corrupt_scandal_sim", 
                          incl.cut=1,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE,
                          dcc = FALSE)
ttout3_sim

#either look manually. No simultaneos subset relations found. Or with formula:
SSR_sim <- intersect(rownames(ttout3_sim$tt) [ttout3_sim$tt$OUT==1],
                                rownames(ttnoout3_sim$tt) [ttnoout3_sim$tt$OUT==1])
SSR_sim
# = 0
# -> no simultaneous subset relations can be found
#4 Incoherent counterfactuals
contcountercrisp2_out3_sim <- findRows(obj=ttout3_sim, type = 2)
contcountercrisp2_out3_sim
# Incoherent rows: 4  7 12 20 23 28 31 32

#to be able to compare them, a new truth table is constructed
ttout3_simenhanced <- truthTable(data=df_set_out2_sim, outcome = "out_removal_sim", conditions = "elite_ide_polar_sim, elite_aff_polar_sim, mass_mobilization_sim, leg_shield_sim, corrupt_scandal_sim", 
                          incl.cut=1,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE,
                          dcc = TRUE)
ttout3_simenhanced
# turn remainder rows to 0
ttout3_simenhanced$tt['4', 'OUT'] <- 0
ttout3_simenhanced$tt['7', 'OUT'] <- 0
ttout3_simenhanced$tt['12', 'OUT'] <- 0
ttout3_simenhanced$tt['20', 'OUT'] <- 0
ttout3_simenhanced$tt['23', 'OUT'] <- 0
ttout3_simenhanced$tt['28', 'OUT'] <- 0
ttout3_simenhanced$tt['31', 'OUT'] <- 0
ttout3_simenhanced$tt['32', 'OUT'] <- 0

#enhanced truth table
ttout3_simenhanced
#write.csv(ttout3_simenhanced$tt, "ttout3_simenhanced.csv")

#enhanced parsimonious solution
psttout3_simenhanced <- minimize(ttout3_simenhanced, include="?", details=TRUE, show.cases=TRUE, row.dom=TRUE, all.sol=TRUE, use.tilde=TRUE)
psttout3_simenhanced 

#enhanced intermediate solution
isttout3_simenhanced <- minimize(ttout3_simenhanced, include="?", dir.exp="-,-,1,0,1", details=TRUE, show.cases=TRUE, row.dom=TRUE,
                         all.sol=TRUE, use.tilde=TRUE)
isttout3_simenhanced
isttout3_simenhanced$i.sol$C1P1$EC 
isttout3_simenhanced$i.sol$C1P1$DC

## Out3: Non-removal from office

ttnoout3_sim <- truthTable(data=df_set_out2_sim, outcome = "~out_removal_sim", conditions = "elite_ide_polar_sim, elite_aff_polar_sim, mass_mobilization_sim, leg_shield_sim, corrupt_scandal_sim", 
                          incl.cut=1,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE,
                          dcc = TRUE)
ttnoout3_sim

isttnoout3_sim <- minimize(ttnoout3_sim, include="?", dir.exp="-,-,0,1,0", details=TRUE, show.cases=TRUE, row.dom=TRUE, 
                         all.sol=TRUE, use.tilde=TRUE)
isttnoout3_sim
isttnoout3_sim$i.sol$C1P1$EC 
isttnoout3_sim$i.sol$C1P1$DC


#Enhanced standard analysis for the absence of the Outcome ~out 

#1 Contradictions to necessity:
# no necessary conditions for Oout -> no contradictions possible
#2 Implausible assumptions
# No implausible assumptions
#3 Simultaneous subset relations 
SSR_sim <- intersect(rownames(ttout3_sim$tt) [ttout3_sim$tt$OUT==1],
                                rownames(ttnoout3_sim$tt) [ttnoout3_sim$tt$OUT==1])
SSR_sim
# = 0
# -> no simultaneous subset relations can be found
#4 Incoherent counterfactuals
contcountercrisp2_nout3_sim <- findRows(obj=ttnoout3_sim, type = 2)
contcountercrisp2_nout3_sim
# Incoherent rows: 4  7 12 20 23 28 31 32
#to be able to compare them, a new truth table is constructed
ttnoout3_simenhanced <- ttnoout3_sim

# turn remainder rows to 0
ttnoout3_simenhanced$tt['4', 'OUT'] <- 0
ttnoout3_simenhanced$tt['7', 'OUT'] <- 0
ttnoout3_simenhanced$tt['12', 'OUT'] <- 0
ttnoout3_simenhanced$tt['20', 'OUT'] <- 0
ttnoout3_simenhanced$tt['23', 'OUT'] <- 0
ttnoout3_simenhanced$tt['28', 'OUT'] <- 0
ttnoout3_simenhanced$tt['31', 'OUT'] <- 0
ttnoout3_simenhanced$tt['32', 'OUT'] <- 0
ttnoout3_simenhanced
#write.csv(ttnoout3_simenhanced$tt, "ttnoout3_simenhanced.csv")

#enhanced parsimonious solution
psttnoout3_simenhanced <- minimize(ttnoout3_simenhanced, include="?", details=TRUE, show.cases=TRUE, row.dom=TRUE, all.sol=TRUE, use.tilde=TRUE)
psttnoout3_simenhanced #enhanced version

#enhanced intermediate solution
isttnoout3_simenhanced <- minimize(ttnoout3_simenhanced, include="?", dir.exp="-,-,0,1,0", details=TRUE, show.cases=TRUE, row.dom=TRUE,
                         all.sol=TRUE, use.tilde=TRUE) 
isttnoout3_simenhanced

```

# Robustness analysis using executive approval (Appendix 14)
```{r}
# load data
data_impeach_approv <- read.csv("data_impeach_approv.csv")

attach(data_impeach_approv)  #attach is used to save a bit time, you do not have to write data_raw$year and so on. 

data_r_approv <- data.frame(code, 
                     out_lower_chamber_vote, 
                     out_removal, 
                     mass_mobilization, 
                     leg_shield, 
                     elite_ide_polar, 
                     elite_aff_polar, 
                     corrupt_scandal, 
                     e_regionpol_7C, 
                     disapproval_not_smoothed, 
                     disapproval_smoothed)

detach(data_impeach_approv)

# --- calibration --- #

df_approv <- data.frame(
  code = data_r_approv$code,
  out_lower_chamber_vote = data_r_approv$out_lower_chamber_vote,
  out_removal = data_r_approv$out_removal,
  elite_ide_polar = calibrate(data_r_approv$elite_ide_polar, 
                              thresholds = 5, 
                              type = "crisp"),
  elite_aff_polar = calibrate(4-data_r_approv$elite_aff_polar, 
                              thresholds = 1.001, 
                              type = "crisp"),
  mass_mobilization = calibrate(data_r_approv$mass_mobilization, 
                                thresholds = 10000, 
                                type = "crisp"),
  leg_shield = calibrate(data_r_approv$leg_shield, 
                         thresholds = 0.1, 
                         type = "crisp"),
  corrupt_scandal = data_r_approv$corrupt_scandal,
  exec_disapprov = calibrate(data_r_approv$disapproval_not_smoothed, 
                               thresholds = 50.0001, 
                               type = "crisp")
)

# subset data
df_set_out1_approv <- df_approv
df_set_out2_approv <- subset(df_approv, out_lower_chamber_vote == 1)

# assign president-years as rownames
df_set_out1_approv <- data.frame(df_set_out1_approv[,-1], row.names=df_set_out1_approv[,1])
df_set_out2_approv <- data.frame(df_set_out2_approv[,-1], row.names=df_set_out2_approv[,1])

# select only relevant variables
df_set_out1_approv_skw <- subset(df_set_out1_approv, select = c(out_lower_chamber_vote, elite_ide_polar, elite_aff_polar, leg_shield, corrupt_scandal, mass_mobilization, exec_disapprov))

df_set_out2_approv_skw <- subset(df_set_out2_approv, select = c(out_removal, elite_ide_polar, elite_aff_polar, leg_shield, corrupt_scandal, mass_mobilization, exec_disapprov))

# check skewness
check_skewness(df_set_out1_approv_skw)
check_skewness(df_set_out2_approv_skw)

# --- Necessity analysis --- #

# Outcome: out_lower_chamber_vote 

cons_out2_approv <- subset(df_set_out1_approv, select = c("elite_ide_polar", "elite_aff_polar", "leg_shield", "corrupt_scandal", "exec_disapprov") )

pof(cons_out2_approv, "out_lower_chamber_vote", df_set_out1_approv, relation = "nec")
pof(1-cons_out2_approv, "out_lower_chamber_vote", df_set_out1_approv, relation = "nec")
pof(cons_out2_approv, "~out_lower_chamber_vote", df_set_out1_approv, relation = "nec")
pof(1-cons_out2_approv, "~out_lower_chamber_vote", df_set_out1_approv, relation = "nec")
#no single necessary condition could be found, neither for the outcome, nor for the non-outcome

#check for SUIN condition
SUIN_out2_approv <- superSubset(df_set_out1_approv, outcome = "out_lower_chamber_vote", conditions = "elite_ide_polar, elite_aff_polar, exec_disapprov, leg_shield, corrupt_scandal", relation = "nec", incl.cut = 0.9, ron.cut = 0.5, cov.cut = 0.6)
SUIN_out2_approv

#absence of outcome
SUIN_noout2_approv <- superSubset(df_set_out1_approv, outcome = "~out_lower_chamber_vote", conditions = "elite_ide_polar, elite_aff_polar, exec_disapprov, leg_shield, corrupt_scandal", relation = "nec", incl.cut = 0.9, ron.cut = 0.5, cov.cut = 0.6)
SUIN_noout2_approv

# Outcome: out_removal

cons_out3_approv <- subset(df_set_out2_approv, select = c("elite_ide_polar", "elite_aff_polar", "exec_disapprov", "leg_shield", "corrupt_scandal"))

pof(cons_out3_approv, "out_removal", df_set_out2_approv, relation = "nec")
pof(1-cons_out3_approv, "out_removal", df_set_out2_approv, relation = "nec")
pof(cons_out3_approv, "~out_removal", df_set_out2_approv, relation = "nec")
pof(1-cons_out3_approv, "~out_removal", df_set_out2_approv, relation = "nec")
#no single necessary condition could be found, neither for the outcome, nor for the non-outcome

#check for SUIN condition
SUIN_out3_approv <- superSubset(df_set_out2_approv, outcome = "out_removal", conditions = "elite_ide_polar, elite_aff_polar, exec_disapprov, leg_shield, corrupt_scandal", relation = "nec", incl.cut = 0.9, ron.cut = 0.5, cov.cut = 0.6)
SUIN_out3_approv

#absence of outcome
SUIN_noout3_approv <- superSubset(df_set_out2_approv, outcome = "~out_removal", conditions = "elite_ide_polar, elite_aff_polar, exec_disapprov, leg_shield, corrupt_scandal", relation = "nec", incl.cut = 0.9, ron.cut = 0.5, cov.cut = 0.6)
SUIN_noout3_approv

# --- Sufficiency analysis --- #

# truth tables for out_lower_chamber_vote
ttout2_approv <- truthTable(data=df_set_out1_approv, outcome = "out_lower_chamber_vote", conditions = "elite_ide_polar, elite_aff_polar, exec_disapprov, leg_shield, corrupt_scandal", 
                          incl.cut=0.8,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE, 
                          dcc = TRUE)
ttout2_approv

# truth tables for out_removal 
ttout3_approv <- truthTable(data=df_set_out2_approv, outcome = "out_removal", conditions = "elite_ide_polar, elite_aff_polar, exec_disapprov, leg_shield, corrupt_scandal", 
                          incl.cut=1,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE,
                          dcc = TRUE)
ttout3_approv


# --- Solutions --- #

## out2: Lower house vote

# intermediate solution isttout2_approv (with directional expectations)
isttout2_approv <- minimize(ttout2_approv, include="?", dir.exp="-,-,1,0,1", details=TRUE, show.cases=TRUE, row.dom=TRUE,
                         all.sol=TRUE, use.tilde=TRUE)
isttout2_approv
isttout2_approv$PIchart
isttout2_approv$i.sol$C1P1$EC #easy counterfactuals
isttout2_approv$i.sol$C1P1$DC


## out3: Removal

# intermediate solution isttout3_approv (with directional expectations)
isttout3_approv <- minimize(ttout3_approv, include="?", dir.exp="-,-,1,0,1", details=TRUE, show.cases=TRUE, row.dom=TRUE,
                         all.sol=TRUE, use.tilde=TRUE)
isttout3_approv
isttout3_approv$PIchart
isttout3_approv$i.sol$C1P1$EC #easy counterfactuals
isttout3_approv$i.sol$C1P1$DC

# --- Enhanced intermediate solution --- #

## Outcome: out_lower_chamber_vote

# Steps:
#1 Contradictions to necessity:
# no necessary conditions for the outcome -> no contradictions possible
#2 Implausible assumptions
# There is no possible implausible assumption in my opinion. 
#3 Simultaneous subset relations 
ttnoout2_approv <- truthTable(data=df_set_out1_approv, outcome = "~out_lower_chamber_vote", conditions = "elite_ide_polar, elite_aff_polar, exec_disapprov, leg_shield, corrupt_scandal", 
                          incl.cut=0.8,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE, 
                          dcc = TRUE)
ttnoout2_approv
# compare if any of the rows that you included into the solution is also part of a possible solution for the non-outcome
ttnoout2_approv
ttout2_approv

#either look manually. No simultaneos subset relations found. Or with formula:
SSR <- intersect(rownames(ttout2_approv$tt) [ttout2_approv$tt$OUT==1],
                                rownames(ttnoout2_approv$tt) [ttnoout2_approv$tt$OUT==1])
SSR
# = 0
# -> no simultaneous subset relations can be found
#4 Incoherent counterfactuals
contcountercrisp2_out2 <- findRows(obj=ttout2_approv, type = 2)
contcountercrisp2_out2
# = 0

#to be able to compare them, a new truth table is constructed
ttout2_approvenhanced <- truthTable(data=df_set_out1_approv, outcome = "out_lower_chamber_vote", conditions = "elite_ide_polar, elite_aff_polar, exec_disapprov, leg_shield, corrupt_scandal", 
                          incl.cut=0.8,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE, 
                          dcc = TRUE)
ttout2_approvenhanced

#enhanced truth table
ttout2_approvenhanced
#write.csv(ttout2_approvenhanced$tt, "ttout2_approvenhanced.csv")

#enhanced parsimonious solution
psttout2_approvenhanced <- minimize(ttout2_approvenhanced, include="?", details=TRUE, show.cases=TRUE, row.dom=TRUE, all.sol=TRUE, use.tilde=TRUE)
psttout2_approvenhanced 

#enhanced intermediate solution
isttout2_approvenhanced <- minimize(ttout2_approvenhanced, include="?", dir.exp="-,-,1,0,1", details=TRUE, show.cases=TRUE, row.dom=TRUE,
                         all.sol=TRUE, use.tilde=TRUE)
isttout2_approvenhanced
isttout2_approvenhanced$i.sol$C1P1$EC 
isttout2_approvenhanced$i.sol$C1P1$DC


## Outcome: out_removal 

# Steps:
#1 Contradictions to necessity:
# no necessary conditions for the outcome -> no contradictions possible
#2 Implausible assumptions
# There is no possible implausible assumption in my opinion. 
#3 Simultaneous subset relations 
ttnoout3_approv <- truthTable(data=df_set_out2_approv, outcome = "~out_removal", conditions = "elite_ide_polar, elite_aff_polar, exec_disapprov, leg_shield, corrupt_scandal", 
                          incl.cut=1,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE,
                          dcc = TRUE)
ttnoout3_approv
# compare if any of the rows that you included into the solution is also part of a possible solution for the non-outcome
ttnoout3_approv
ttout3_approv

#either look manually. No simultaneos subset relations found. Or with formula:
SSR <- intersect(rownames(ttout3_approv$tt) [ttout3_approv$tt$OUT==1],
                                rownames(ttnoout3_approv$tt) [ttnoout3_approv$tt$OUT==1])
SSR
# = 0
# -> no simultaneous subset relations can be found
#4 Incoherent counterfactuals
contcountercrisp2_out3 <- findRows(obj=ttout3_approv, type = 2)
contcountercrisp2_out3
# Incoherent rows: 3  7 19 23

#to be able to compare them, a new truth table is constructed
ttout3_approvenhanced <- truthTable(data=df_set_out2_approv, outcome = "out_removal", conditions = "elite_ide_polar, elite_aff_polar, exec_disapprov, leg_shield, corrupt_scandal", 
                          incl.cut=1,
                          n.cut=1, 
                          sort.by="incl, n", 
                          complete=TRUE, 
                          show.cases=TRUE,
                          dcc = TRUE)
ttout3_approvenhanced
# turn remainder rows to 0
ttout3_approvenhanced$tt['3', 'OUT'] <- 0
ttout3_approvenhanced$tt['7', 'OUT'] <- 0
ttout3_approvenhanced$tt['19', 'OUT'] <- 0
ttout3_approvenhanced$tt['23', 'OUT'] <- 0

#enhanced truth table
ttout3_approvenhanced
#write.csv(ttout3_approvenhanced$tt, "ttout3_approvenhanced.csv")

#enhanced parsimonious solution
psttout3_approvenhanced <- minimize(ttout3_approvenhanced, include="?", details=TRUE, show.cases=TRUE, row.dom=TRUE, all.sol=TRUE, use.tilde=TRUE)
psttout3_approvenhanced 

#enhanced intermediate solution
isttout3_approvenhanced <- minimize(ttout3_approvenhanced, include="?", dir.exp="-,-,1,0,1", details=TRUE, show.cases=TRUE, row.dom=TRUE,
                         all.sol=TRUE, use.tilde=TRUE)
isttout3_approvenhanced
isttout3_approvenhanced$i.sol$C1P1$EC 
isttout3_approvenhanced$i.sol$C1P1$DC


```

# Regression analysis (Appendix 15)
```{r}

library(stargazer)

# load data 
data_raw <- read.csv("data_impeach_main.csv")

data <- data_raw %>%
  mutate(
    elite_aff_polar_rev = -elite_aff_polar,  # reverse coding
    log_mass_mobilization = log1p(mass_mobilization)  # log(1 + x)
  )

# Outcome 1: Lower house vote 

model1 <- glm(out_lower_chamber_vote ~ elite_ide_polar + elite_aff_polar_rev +
                corrupt_scandal + leg_shield + 
                log_mass_mobilization, data = data, family = binomial)

summary(model1)

model2 <- glm(out_lower_chamber_vote ~ elite_ide_polar * corrupt_scandal + elite_aff_polar_rev +
                + leg_shield + 
                log_mass_mobilization, data = data, family = binomial)

summary(model2)

model3 <- glm(out_lower_chamber_vote ~ elite_ide_polar + elite_aff_polar_rev * corrupt_scandal +
                   log_mass_mobilization  + leg_shield  
                , data = data, family = binomial)

summary(model3)

model4 <- glm(out_lower_chamber_vote ~ elite_ide_polar * log_mass_mobilization + elite_aff_polar_rev + 
                   corrupt_scandal + leg_shield  
                , data = data, family = binomial)

summary(model4)

model5 <- glm(out_lower_chamber_vote ~ elite_ide_polar + elite_aff_polar_rev * log_mass_mobilization + 
                   corrupt_scandal + leg_shield  
                , data = data, family = binomial)

summary(model5)

# Outcome 2: Removal from office

model6 <- glm(out_removal ~ elite_ide_polar + elite_aff_polar_rev +
                corrupt_scandal + leg_shield + 
                log_mass_mobilization, data = data, family = binomial)

summary(model6)

model7 <- glm(out_removal ~ elite_ide_polar * corrupt_scandal + elite_aff_polar_rev +
                + leg_shield + 
                log_mass_mobilization, data = data, family = binomial)

summary(model7)

model8 <- glm(out_removal ~ elite_ide_polar + elite_aff_polar_rev * corrupt_scandal +
                   log_mass_mobilization  + leg_shield  
                , data = data, family = binomial)

summary(model8)

model9 <- glm(out_removal ~ elite_ide_polar * log_mass_mobilization + elite_aff_polar_rev + 
                   corrupt_scandal + leg_shield  
                , data = data, family = binomial)

summary(model9)

model10 <- glm(out_removal ~ elite_ide_polar + elite_aff_polar_rev * log_mass_mobilization + 
                   corrupt_scandal + leg_shield  
                , data = data, family = binomial)

summary(model10)

# Output results

m_list <- list(model1, model2, model3, model4, model5, model6, model7, model8, model9,model10)

stargazer(m_list, title = "Impeachment Hurdle 1 (models 1 - 5) & Hurdle 2 (models 6 - 10)", type = "text", model.numbers = TRUE)

```

